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Abstract 

I consider the appearance of shocks in hyperbolic formalisms of General Rel- 
ativity. I study the particular case of the Bona-Masso formalism with zero 
shift vector and show how shocks associated with two families of characteristic 



. fields can develop. These shocks do not represent discontinuities in the geom- 



etry of spacetime, but rather regions where the coordinate system becomes 
pathological. For this reason I call them 'coordinate shocks'. I show how 
one family of shocks can be eliminated by restricting the Bona-Masso slicing 
condition dta = — a;^/(a)tri^ to the case / = 1 + kjo? , with k an arbitrary 



^ . constant. The other family of shocks can not be eliminated even in the case 

of harmonic slicing (/ = 1). I also show the results of the numerical evolution 
of non-trivial initial slices in the special cases of a flat two-dimensional space- 
time, a flat four-dimensional spacetime with a spherically symmetric slicing, 
and a spherically symmetric black hole spacetime. In all three cases coordi- 
nate shocks readily develop, confirming the predictions of the mathematical 
analysis. Although I concentrate in the Bona-Masso formalism, the phenom- 
ena of coordinate shocks should arise in any other hyperbolic formalism. In 
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particular, since the appearance of the shocks is determined by the choice of 

gauge, the results presented here imply that in any formalism the use of a 

harmonic slicing can generate shocks. 
04.20.Ex,04.20.Dm 
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I. INTRODUCTION 



In the last few years there has been a renewed interest in the study of initial- value formu- 
lations of General Relativity This interest has been motivated mainly by the desire of 
rewriting the Einstein system of evolution equations in an explicit hyperbolic form, so that 
it can be solved numerically using modern high resolution methods from fluid dynamics |]TU 



One can separate the new hyperbolic formalisms according to the way in which they 
treat the evolution of the lapse function a. Some formulations assume the existence of an 
arbitrarily prescribed gauge, i.e. the lapse is an arbitrary function of spacetime known a 
'priori [H,^ ('prescribed gauge' formalisms). Other formulations include the lapse function 
as part of the system of dynamical variables, and postulate for it an evolution equation that 
guarantees the hyperbolicity of the whole system, geometry plus gauge ('hyperbolic 
gauge' formalisms). The resulting formalisms remain hyperbolic for any prescribed shift 
vector. 

Prescribed gauge formalisms, though certainly useful theoretically, might have a limited 
applicability in Numerical Relativity simply because there is no recipe that can give us the 
a priori form of the lapse except in trivial cases (for example a = 1). Hyperbolic gauge 
formalisms on the other hand, by allowing the lapse functon to adapt itself to the evolution 
of the geometry while maintaining the hyperbolic structure of the system of equations, would 
appear to be much more promising. 

Hyperbolic gauge formalisms, however, are probably more susceptible to a problem that 
seems to have been overlooked until know. By rewriting the whole evolution system (gauge 
plus geometry) in hyperbolic form, they open up the possibility of running into a well known 
non-linear effect associated with hyperbolic systems: the appearance of shocks. Here I use 



^To date, in all these hyperbolic formalisms the shift vector is assumed to be known a •priori. The 
author is aware, however, of some efforts to find evolution equations for the shift that will keep the 
whole system hyperbolic [O]. 
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the term 'shock' in a somewhat loose form to refer to a discontinuous solution that develops 
from smooth initial data, without worrying about the existence of weak solutions or jump 
conditions. 

The fact that in vacuum General Relativity one can have shock fronts is well 
known ||T2|-|r^. By shocks fronts, however, one generally understands discontinuities in 



the curvature of spacetime present in the initial data that propagate with the speed of light. 
In the theory of non-linear hyperbolic equations such solutions are not considered proper 
shocks, but are called instead 'contact discontinuities'. Here, however, I will consider the 
existence of discontinuous solutions that arise from smooth initial data even in a flat space- 
time. Clearly those solutions do not correspond to a physical discontinuity in the geometry 
of spacetime. Instead the discontinuities indicate regions where our coordinate system be- 
comes pathological: the time slices can become non-smooth, or a spatial coordinate might 
map a flnite proper distance to an infinitesimal interval. It is for this reason that I shall 
refer to them as 'coordinate shocks'. 

Even though modern high resolution numerical methods can deal with the presence of 
shock waves, clearly the appearance of coordinate shocks is something that must be avoided. 
In the first place, coordinate shocks create completely artificial discontinuities in solutions 
that otherwise represent perfectly smooth geometries. Not only that, but since in general our 
gauge conditions are not obtained from a conservation law, we will not have an analogue of 
the 'weak solutions' to such laws. This means that after a shock forms our gauge conditions 
will just break down, and even if the numerical solution remains well behaved, it will not 
have any clear physical meaning. In particular, as the numerical mesh is refined, the solution 
will not converge after the formation of the shock. 

In this paper I will concentrate in one particular hyperbolic formalism of General Rela- 
tivity, the Bona-Masso (BM) formalism and I will show how these coordinate shocks 
can and do indeed develop even in very simple situations. 
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II. THE BONA-MASSO FORMALISM 



In this section I will make a brief introduction to the BM hyperbolic formalism for 
General Relativity. I will use the most recent form of this formalism as presented in [Q. 

Let us start from the standard 3+1 formulation of General Relativity of Arnowitt, Deser 
and Misner (ADM) [0,2C]. The evolution equations for the metric Qij and extrinsic curva- 
ture Kij are: 



{dt - Cp) Qij = -2 a Kij , 



(la) 

(lb) 



where a is the lapse function, (3^ the shift vector, and where Rf^ and i?-^"* represent 
the components of the Ricci tensor for the spatial hyper- surfaces and for the full spacetime 
respectively. In what follows I will restrict myself to the case of zero shift vector. The ADM 
equations then reduce to: 



dt 9ij = -2 a Kii 



dtKij = -ViVja + a + tiK Kij - 2K,kK^ - R 



(2a) 
(2b) 



In order to obtain a system that is first order in space we introduce the following quan- 
tities: 



Au = dklna 



(3) 



The evolution equation for Kij can then be rewritten as: 



dt Kij + dk (a X^ij^ = aS, 



(4) 



where we have defined: 



D^, + ^ 51 {A, +2Vj- D,„r) + ^ 5|(A, + 2 - A„ 



(5) 



with: 



Vk = Dkni^ — D'^mk ■ (6) 

The source term Sij in equation involves only the fields themselves and not their 
derivatives: 

+ [A' - 2Dj^) (A,fc + D,,k) + A, [Vj - \ D,k') + A, [v, - \ A/) • (7) 

We also need an evolution equation for the Dkij- This we obtain by taking the spatial 
derivative of the evolution equation for Qif 

^tDk^J + dk{aKi,) = . (8) 

The quantities Vk defined in are very important. Their evolution equation can be 
obtained from (^). In order to ensure hyperbolicity, however, it is crucial to modify the 
resulting equation using the momentum constraints to obtain: 

dtVk = aPk , (9) 

where: 

Pk = Gl + A^ [K^ - 5™trir) + (A,n" - 5fc /^^a") 

- 2K^^ {d"^\ - 51 D^^"^) , (10) 

and where G^^^ is the Einstein tensor of the spacetime. 

The quantities Vk are now considered independent, and equation (||) becomes an alge- 
braic constraint that must be satisfied by the physical solutions. 

Finally, we need evolution equations for the lapse a and its derivative Ak , i. e. we need 
to choose a slicing condition. In the BM formalism the following slicing condition is used: 

dta = -a^ f{a) tiK , (11) 

with /(a) > but otherwise arbitrary. 
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The complete system of evolution equations then takes the form: 

dta = -a'^ f (a) trK , (12a) 

dtgij = -2aKi, , (12b) 

and: 

dtAk+ dk{aftiK) =0 , (13a) 

^tDk^, + dk{aKij) =0 , (13b) 

dtKij + dk (aA^-'jj) = a Sij , (13c) 

dtVk = aPk . (13d) 

To study the characteristic structure of the system of equations (|T^) and (|TB|) we choose 
a fixed space direction x and consider only derivatives along that direction. It can then be 
shown that the system is hyperbolic with the following structure: | 

• 25 fields propagate along the time lines (zero speed). These fields are: 

V„ Ax - fD^m} ix' ^ x) . (14) 

• 10 fields propagate along the physical hght cones with speeds: 

A'± = ±aJg^ . (15) 



These fields are: 

.± = K,x' ± ip^ix' + 5f KV^?^^) {x' + x) . (16) 



^Here I use the term hyperbolic in the weak sense to mean that the characteristic matrix of the 
system has real eigenvalues. It should be noticed that this weak form of hyperbolicity does not 
guarantee that the system can be diagonalised. A crucial feature of the BM formalism is that even 
though it is only weakly hyperbolic, it can in fact be diagonalised as long as / > . 
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2 fields propagate with the 'gauge speeds': 

A^± = ±a^[Jg- 



(17) 



They are: 

= ^tiK ± (A^ + 2\/7^^^) . (18) 



III. MATHEMATICAL ANALYSIS OF THE NON-LINEARITIES 

Here I will try to understand the nonlinear it ies present in our system of equations, trying 
in particular to determine whether shocks can develop. From the discussion in the previous 
section it is clear that the system of evolution equations in the BM formalism has the 
following structure: 

dtUi = Pi i e {1, ...,N^} , (19a) 

dtv^ + d^Fi = te {l,...,iV4 . (19b) 

The fluxes Fi that appear in the above equations have the form: 

i=i 

where the coefficients Mij are functions of the m's but not of the f 's. 

Let us now call Xi the eigenvalues and Gj the corresponding eigenvectors of the Jacobian 
matrix Mij = d Fi / dvj . Let us also introduce the matrix R = [ei|e2| • • ■ |eAr^] of column 
eigenvectors. The eigenfields Wi are then defined by: 

Y = Rw w = R-^v . (21) 



A given eigenfield Wi is called 'linearly degenerate' [21| if the following condition holds: 



|^.£§^|l^.V„A..e.^O. (22) 

owi p{ dvj dwi 



Since in our case the A's don't depend on the w's, it is obvious that all the eigenfields 
are linearly degenerate. 

In the case of systems of conservation laws where the sources vanish, linear degeneracy 
is enough to guarantee that no shocks will form. However, when the sources are non-zero, 
this is not true anymore. This is easy to see if we consider for a moment the prototype of 
non-linear hyperbolic equations. Burgers' equation: 

dtu + udxU = . (23) 

If we now define: 

V := d^u , (24) 
then we can rewrite equation (|23| ) as the system: 



dtU = —uv , (25a) 
dtv + dx (uv) = . (25b) 

This has precisely the form ([T9|) . The only eigenvalue turns out to be equal to u which is 
clearly independent of v. By the definition above the system is linearly degenerate. However, 
it is clearly non-linear and will generate shocks since it is only Burgues' equation in disguise. 
The nonlinearities have now been buried in the sources. 

Clearly the condition that must be imposed to guarantee that no shocks will develop is 
that a given eigenvalue Aj should not be affected by changes in the corresponding eigenfield 
Wi. The condition for linear degeneracy ( p2D asks for the eigenvalue not to be explicitly 
dependent on its associated eigenfield. In the presence of sources, however, the coupling 
can introduce an indirect dependency. In order to study this dependency let us consider the 
time evolution of Aj : 

\i = dtK = dt Uj = V„ Ai ■ p . (26) 

Now, we want this time derivative to be independent of the eigenfield Wi. 
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dXi d 



(V, Ai ■ p) = . 



(27) 



dwi dwi 

I shall call this condition 'indirect linear degeneracy' and I will refer to condition ( P^ as 
'explicit' or 'direct' linear degeneracy. 

If we assume that the condition for explicit linear degeneracy holds, then the condition 
for indirect linear degeneracy can be reduced to: 



V„, A,; ■ = V^. A, ■ 2^ TT— -TT^ = 



dWi 



r{ dvj dwi 



(28) 



which can be rewritten as: 



V«Ai ■ (ej ■ V„) p = 



(29) 



This condition must supplement the condition for explicit linear degeneracy ( p2|) if we 
want to guarantee that no shocks will develop. 



Let us now apply the previous condition to the BM system of evolution equations (|T2D 
and ([l3|) . From the discussion of the previous section it is clear that, on a given spatial 
direction x we only have the following non-trivial eigenvalues: 



A'± = ±aJg-- , Xf± = ±aJf 



The time derivative of A^± will then be: 
A'± = ±A^ 



- 9ta + dt g 
a 2g^^ 

gXm gXn 

a 1 g^^ 



Using now equations ([T2|) we find: 



A'± = ±aA'± [K''''lg'''' - ftiK 



(30) 



(31) 



(32) 



Now, from the definitions of wi and Wf (equations ([16D and ([T8|)) we can easily find that: 



trK 



2v7 



(33) 
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and (j9, q ^ x): 



^ gXX ^ (^gXp gXg _ 



v7 



pq- 



(34) 



Substituting these results back in the expression for A'± we find: 



^ (1 - /) {u.' 



In the same way we find for the time derivative of A-^^: 



(35) 



A/± = ±aXL 



j^xx^gxx _ ^ ^ tiK 



(36) 



where /' = f ■ Substituting again the expression for K^^ and tri^ in terms of the 
eigenfields we find: 



A/+ = ± 



(l - / - « f/2) {w^^ + wf_) 



/ 

pq- 



(37) 



Equations ( ^5|) and (|3^) are very important resuhs. Consider first the situation for X^^^. 
If we want A-^ ± to be independent of w±. , and hence satisfy the condition for indirect linearly 
degeneracy, we must clearly ask for: 



1 - / - a/72 = 



This differential equation can be easily solved to give: 



/(«) = ! + k/a'^ 



(38) 



(39) 



with k an arbitrary constant. We must in fact take > in order to ensure that we will 
have / > for all a > . 

We have then show that the function / must have the form in order to guarantee 
that the eigenfields w£ will not generate shocks. Notice that if we take k = the condition 
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reduces to that of harmonic shcing, i.e. for harmonic shcing the eigenfields w± do not 
generate shocks. 

Consider now the situation for A'±. From equation 
to be independent of w^^^ we must have: 

This condition is very restrictive. In particular, it is impossible to satisfy with a diagonal 
metric. We then reach the conclusion that in the general case, the eigenfields Wgp-t can 
always generate shocks. Notice how this result is independent of the value of /, it will 
therefore remain true even in the case of harmonic slicing. 

One must stress here the fact that we haven't actually proved that shocks will indeed 
develop. Whether they do or not in any particular case should depend in a critical way on 
the form of the initial data. 

In the following sections I will consider some examples that show how coordinate shocks 
can indeed develop even in very simple cases. 



35) it is clear that if we want A'- 



IV. FLAT TWO-DIMENSIONAL SPACETIME 

A. Evolution equations 

As a first example, consider a flat two-dimensional spacetime (a '1+1' spacetime) with co- 
ordinates {t, x} . Notice that these coordinates do not have to correspond to the Minkowski 
coordinates {xmj^m} ? so we can have a non-trivial evolution even though the spacetime is 
flat. Since we only have one spatial dimension, I will simplify the notation in the following 
way: 

g := g^^ , A := , D := D^^^ , K := K^^ . (41) 
Notice that the variable is identically zero. 
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The system of evolution equations (|T^) and (|T3D reduces in this case to: 

dta = -a^fK/g , (42a) 
dtg = ~2aK , (42b) 



and: 



dtA + d,{afK/g) =0 , (43a) 

dtD + d,{aK) =0 , (43b) 

dtK + [a a) = a/g [aD - K^) . (43c) 

The characteristic structure of this system is very simple: 

• There are 3 fields that propagate along the time lines (speed zero). These fields are: 

{a,g,A-fD/g} . (44) 

• The 2 remaining fields propagate with the 'gauge speeds': 



A^± = ±a^f/g . (45) 

They are: 



w^± = K/g ± A/^g . (46) 

Notice how there are no fields propagating along the physical light cones. According to 
the discussion of the previous section, we should then expect shocks only when condition (|39D 
is violated. 

B. Numerical simulations 

Since we are dealing with a flat spacetime, the only way to obtain a non-trivial evolution 
is to start with a non-trivial initial slice. I will therefore consider an initial slice given in 
terms of Minkowski coordinates {xm, thi} as: 

13 



tM = h{xM) ■ (47) 

I will assume that the dynamical spatial coordinate x coincides initially with the 
Minkowski spatial coordinate xm- It is then not difficult to show that the initial metric 
g and extrinsic curvature K are given by: 

g = l - h'^ , (48a) 
K = -h"/^. (48b) 

The initial value of D can be obtained directly from its definition in terms of g. The 
initial lapse is taken to be equal to 1 everywhere, which implies that A = 0. 
In all the simulations shown here, the function h{x) has a Gaussian profile: 

h{x) = H expj-^^^^l , (49) 

with {H, a, Xc} constants. The particular values of {H, a} used in the simulations presented 
here are: 

H = 5 , a = 10 . (50) 

I have also always taken the initial perturbation to be centered around Xc = 150. The 
initial values of all the variables can be seen in Figure [l|. All the results presented below 
where obtained using a time step of At = 0.125 and a spatial increment of Ax = 0.25. 

In all the simulations, the evolution proceeds at ffist in a similar way: The initial per- 
turbation in g , D and K gives rise to perturbations in a and A . These perturbations 
rapidly develop into two separate pulses traveling in opposite directions with a speed ~ . 
What happens later depends crucially on the form of the function f{a). 

For harmonic slicing (/ = 1), the pulses remain smooth as they move away. Once the 
pulses are gone, the lapse, the metric, and the variables A and D return to their initial 
values, and the extrinsic curvature becomes 0. Figure ^ shows the values of the variables at 
t = 100. 
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When / is a constant larger than 1, the pulses do not remain smooth and shocks develop. 
In fact, we have two shocks developing in each pulse, one in front of it and one behind it. At 
those points, the lapse and the metric develop large gradients, while the extrinsic curvature 
and the variables A and D develop very tall and narrow spikes. Figure |^ shows the values 
of the variables at t = 75 in the particular case when / = 1.69. 

When / is a constant smaller than 1, a single shock develops in the middle of each pulse. 
Appart from this, the situation is very similar to the case / > 1. Figure § shows the values 
of the variables at t = 75 in the particular case when / = 0.49. 

Finally, when / is of the form (|39|) , no shocks develop in agreement with the predictions. 
The pulses remain smooth and move away with a speed ~ Vl + k. Figure || shows the values 
of the variables at t = 70 in the particular case when / = 1 + 1 / . 

As a final comment, it should be mentioned that I have performed similar simulations 
with many different values of the amplitude H and width a of the gaussian profile of the 
initial slice. When / is not of the form (|39| ) shocks apparently always develop, though at 
different times. The crucial feature that seems to determine the time of shock formation is 
the maximum absolute value of the extrinsic curvature K. For small values of K , shocks 
take a long time to appear, whereas for large values they develop very rapidly. 

V. SPHERICALLY SYMMETRIC VACUUM SPACETIME 

A. Evolution equations 

As a second example, consider a spherically symmetric four-dimensional vacuum space- 
time. Let us introduce the coordinate system {t, r, 9, 0} . The only independent dynamical 
variables will then be: 

|q;, g^r, gee, ^r, D^rr, Dree, K„, Kqo, V^j . (51) 
The system of evolution equations ([1^) and (pTSD reduces now to: 
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and: 



with: 



and: 



dta — — f trK , 

dt Qrr = -'^a Krr , 

dtgee^-'^oi Keg , 



dtAr + dr {aftxK 

dt Drrr + dr (a Krr 

dt Dree + dr (a Kee 

dtKrr + dr (a X^.^. 

dtKee + dr {oiXee 



= , 
= , 
= , 

= a See , 



dfVr = a Pr . 



A;, = A + 2Vr- 2 Dree/ gee , 
^ee — Dree/ 9rr ■, 



Srr = Krr Kge/ Qee — Krr/ 9rr) + A {Drrr/ Qrr " 2 Dree/ 9ee) 

+ 2 Dree/gee i^Drrr/grr — Dree/gee^ +2ArVr , 

See — Krr Kee/ grr — Drrr Dree/ grr^ + 1 , 



ArKee - Dree [Kee/gee - Krr/gr^ 



(52a) 
(52b) 
(52c) 



(53a) 
(53b) 
(53c) 
(53d) 
(53e) 
(53f) 



(54a) 
(54b) 



(55a) 
(55b) 

(55c) 
(55d) 



We also have the following algebraic constraint that must be satisfied by the physical 
solutions: 



Vr — Dree /gee ■ 



The characteristic structure of this system turns out to be: 



(56) 
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• 5 fields propagate along the time lines (speed zero). These fields are: 

{a, 9rr, gee, K, A - f Dr'^m} ■ (57) 

• 2 fields propagate along the physical light cones with speeds: 

A'± = ±a/^r . (58) 

These fields are: 

ty'± = Vg^ Kee ± Dree ■ (59) 

• 2 fields propagate with the 'gauge speeds': 

A^± = ±a//7^. (60) 

They are: 

w^± = \fTg^r trX ± (A^ + 2 Vr) . (61) 

Notice how we now have both fields propagating with the speed of light and fields prop- 
agating with the gauge speed. We should then expect to see two different types of shocks 
forming. In particular, shocks produced by the w^± fields can be expected always, even for 
harmonic slicing. 

B. Numerical simulations for a flat spacetime 

Again, since we are dealing with flat spacetime, the only way to obtain a non-trivial 
evolution is to start with a non-trivial initial slice. I will therefore consider an initial slice 
given in terms of Minkowski coordinates {r^, tu} as: 

tM ^ h {tm) . (62) 
17 



I will assume that the dynamical radial coordinate r coincides initially with the 
Minkowski radial coordinate tm- It is then not difficult to show that the initial metric 
{grr,g9e} and extrinsic curvature {KrrjKee} are given by: 



grr = I — h' , (63a) 

gee = r'^ , (63b) 

Krr = — h" / y/g^ , (63c) 

Kee = - rh' / y/g^ . (63d) 



The initial values of {Drrr, -Dree, K} can be obtained directly from their definitions in 
terms of the metric. The initial lapse is taken to be equal to 1 everywhere, which implies 
that Ar = 

In all the simulations shown here, the function h{r) has a Gaussian profile: 

h{r) = H exp|-^^^-^| , (64) 

with {if, 0", Tc} constants. The particular values of {H, a} used in the simulations presented 
here are: 

H = Ih , (J = 20 , (65) 

and I have taken the initial perturbation to be centered around = 300. The initial values 
of all the variables can be seen in Figure ^ The results presented below where obtained 
using a time step of At = 0.1 and a spatial increment of Ax = 0.2. 

In all the simulations the evolution proceeds at first in a similar way: the initial pertur- 
bations in {grr, Drrr, Krr, Kgg} give risc to perturbatious in {a, ggg, Ar, Dree, K}- These 
perturbations develop into separate pulses traveling in opposite directions with a speed 
~ vT- The pulses are not symmetric any more since clearly the in going and out going 
directions are not equivalent. 

Consider first the case of / > 1. As the evolution proceeds, shocks develop in both 
pulses. These shocks are similar to those found in the 1+1 case: two shocks develop in each 
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pulse, one in front of it and one behind it. At those points {a, grr, Dree, K00, Vr} develop 
large gradients, while {Ar, Drrr, K„} develop tall and narrow spikes. The angular metric 
component gee in contrast develops sharp corners. Figure |^ shows the values of the variables 
at t = 70 in the particular case when / = 1.69. 

When / < 1, we again find results that are similar to the 1+1 case: a single shock 
develops in each pulse. Again, at the shock {a, g^, Dj-ee, K00, Vr} develop large gradients, 
{Ar, Drrr, Krr} dcvclop spikcs and gee develops sharp corners. Figure || shows the values 
of the variables at t = 70 in the particular case when / = 0.49. 

The most interesting case is that of harmonic slicing (/ = 1). In contrast to the 1+1 
case, shocks still develop here. The shocks, however, have a different structure indicative 
of their different origin: the variables {Ar, Drrr, Krr} now develop large gradients, while 
{a, grr, Dree, Kee, K-} develop sharp spikes. The angular metric component gee also seems 
to develop a large gradient, thought this gradient is less sharp than that found in other vari- 
ables. This is easy to understand geometrically: any discontinuity in gee must necessarily 
be accompanied by an infinite value of grr (we must jump a finite radial distance in an 
infinitesimal interval). The shocks are clearly visible in the in going pulse, but don't seem 
to be present in the out going pulse. Figure |^ shows the values of the variables at t = 70 
for harmonic slicing. 

Again, I have performed similar simulations for different amplitudes and widths of the 
gaussian profile of the initial slice. The shocks are always there. The only exception seems 
to be that for very small amplitudes of the initial perturbation the pulses moving in the in 
going direction may not have enough time to develop before they reach the origin. Also, 
for the case of harmonic slicing, shocks in the out going direction never seem to form, no 
matter how large the initial perturbation might be. 
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C. Numerical simulations for a black hole spacetime 



In all the previous examples I have restricted myself to a flat spacetime. Since this is a 
very special case one might think that the shocks that we have found are just an artifact of 
the flatness. To show that this is not the case, I will now consider a spherically symmetric 
black hole spacetime. 

To flnd adequate initial data I start from a Schwarzschild slice with spatial metric: 

where is the Schwarzschild radial coordinate and c?f2^ = dO'^ + sin^^ dcj)^. 

In order to eliminate the singularity at = 2M, I will deflne a new radial coordinate r 
that measures proper distance along the slice. The coordinates and r will be related by: 

Tg + Tjir) 



r = 7] (vs) + M In 
with: 



r]{rs) = {rl - 2Mr,)'^' . (68) 

Notice that even though (^) can not be inverted analytically to flnd rs(r), it can easily 
be inverted numerically to arbitrarily high accuracy. 
The new metric will now have the form: 

= ^ (r,(r))^ dQ^ . (69) 

It is easy to see that the Schwarzschild slice has zero extrinsic curvature, so our initial 
data will be: 

grr = l , (70a) 

gee = r', , (70b) 

Krr = , (70c) 

Kee = . (70d) 
20 



(67) 



Now, if we use this initial data directly we will not see any shocks develop. This is 
known since the BM formalism has been used before to solve this problem and no shocks 
have been observed [0]. The reason why shocks don't develop is that they are a consequence 
of transport and as such they should only develop when we have wave propagation, either 
in the form of real gravitational waves, or in the form of pure gauge waves. The static black 
hole problem has no gravitational waves, and the initial data given above does not give rise 
to gauge waves either. 

In order to introduce gauge waves into our problem, I will consider an initial slice given 
in terms of Schwarzschild time tg in the following way: 

ts = h{r) . (71) 

It is not difficult to show that the new slice will have the following metric components: 

grr = l- [as h'^ , (72a) 
9ee = rl , (72b) 

where as is the Schwarzschild lapse function: 

,1/2 



a. 



(l - 2M/r,) . (73) 
The components of the extrinsic curvature for this slice can now be shown to be: 



a.. 



h" + a'^h' (2 - {ash'f) , (74a) 



Kee = -airs h' / ■ (74b) 

As before, the initial values of {Drrr^ Dree, K} can be obtained directly from their defi- 
nitions in terms of the metric. The initial lapse is again taken to be equal to 1 everywhere, 
which implies that = 0. 

For the function h{r) I will again use a Gaussian: 

h{r) = H expl-^^^^^] , (75) 



a" 
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with {H, a, Tc} constants. 

In order to see the development of the shocks clearly, I will consider simulations where 
the center of our perturbation Tc is out in the wave zone. 

All the simulations I have carried out proceed in a similar way. At the throat of the 
wormhole we find what we expect for a black hole spacetime: the lapse collapses and the 
metric component Qrr grows rapidly. Out in the wave zone, the disturbance behaves in 
the same way as it did in flat spacetime: the initial perturbations in {grr, Drrr, K„, Kqq} 
give rise to perturbations in {a, ggo, Ar, Dree, K}, these then develop into separate pulses 
traveling in opposite directions with a speed ~ v7 . 

In all cases, the traveling pulses develop shocks that have very similar characteristics to 
those that we found in the fiat case. Here I will only show the results found in the case of 
harmonic slicing / = 1. The particular values of {H,a} used in this simulation are: 

H = 5 , a = 5 . (76) 

I have also taken the initial perturbation to be centered around Tc = 50, and the mass of 
the black hole to be M = 1. The results presented here where obtained using a time step 
of At = 0.025 and a spatial increment of Ax = 0.05. The initial values of all the variables 
can be seen in Figure 



Figure [III shows the values of the variables at t = 15. Notice how around the throat the 
lapse and the angular metric component gee have collapsed, while the radial metric compo- 
nent grr has grown to a very large value. The interesting region for our purposes, however, 
is away from the throat. We can clearly see the two pulses resulting from our initial per- 
turbation. The pulse moving inwards has developed a shock: the variables {Ar, Drrr, Krr} 
have developed large gradients, while {a, grr. Dree, Kee, K} have developed sharp spikes. 
The angular metric component gee has also developed a large gradient. 
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VI. DISCUSSION 



I have introduced a general approach to the study of shock development in hyperbolic 
systems of equations with sources. I have shown that the usual condition of explicit linear 
degeneracy (direct linear degeneracy) must be supplemented with a new condition which I 
have called 'indirect linear degeneracy' in order to guarantee that no shocks will develop. 

I have applied this condition of indirect linear degeneracy to the BM hyperbolic formalism 
of General Relativity in the case of a zero shift vector. My analysis has shown how two 
distinct families of characteristic fields can give rise to shocks. Numerical simulations have 
confirmed these predictions in the simple cases of a flat two-dimensional spacetime, a flat 
four-dimensional spacetime with spherically symmetric slices, and a spherically symmetric 
black hole spacetime. 

The appearance of shocks that develop from smooth initial data in vacuum General 
Relativity comes as a great surprise. These shocks, however, do not represent discontinuities 
in the geometry of spacetime, but indicate instead regions where our coordinate system 
becomes pathological. It is for this reason that I refer to them as 'coordinate shocks'. 

Of the two families of coordinate shocks found, one can be completely eliminated by 
choosing a BM gauge function f{a) of the form: 



with A; > an arbitrary constant. For A; > 0, however, this form of the function / will not 
be very useful in spacetimes with large curvatures. The reason for this is easy to see. Even 
thought the condition will prevent the formation of shocks, it implies an evolution equation 
for the lapse of the form: 



Clearly, in a region where the lapse has collapsed to a very small value we will have: 



/ (a) = 1 + A; / , 



(77) 




(78) 



dta ~ — k trK . 



(79) 
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If tri^ > 0, there is nothing to prevent the lapse from becoming negative (this can in fact 
happen very easily in black hole simulations). We are then led to the conclusion that the 
only value of / that will prevent the first family of shocks from developing without carrying 
the risk of leading to a negative lapse is f = 1, i.e harmonic slicing. 

The second family of shocks, on the other hand, is independent of the form of / and 
arises even for harmonic slicing. This is a very unexpected result. After all, this is precisely 
the slicing used to prove the theorems of existence and uniqueness of solutions in General 
Relativity |ll-|l8|. Since at a shock the differential equations break down, one would expect 



the theorems to forbid such solutions. We must remember, however, that these theorems 
can only be proved locally, they can not therefore rule out a shock that develops after a 
finite time. 

It must be stressed that the violation of indirect linear degeneracy is not a sufficient 
condition for the development of shocks. The choice of initial data will have a crucial effect 
in whether or not shocks actually develop. In particular, since shocks are a consequence of 
transport, they should only develop when we have wave propagation, either in the form of 
real gravitational waves, or in the form of pure gauge waves as was shown in the examples 
presented here. Of course, in the simple cases considered in this paper one can easily find 
initial data that does not produce shocks: for the fiat spacetimes one can just take a fiat 
initial slice, while for a black hole spacetime we can start from an unperturbed Schwarzschild 
slice. In the more general case, however, it might be difficult to find such benign initial data, 
or even to prove that it exists at all. 

One more important point should be made here. Since the shocks that I have found 
arise in the case of a zero shift vector, they must necessarily indicate a break down of 
the slicing condition. That is, the shocks represent places where the spatial hyper- surfaces 
become non-smooth. Since the presence of a shift vector can not alter the geometry of these 
hyper-surfaces, the shocks found here must appear for any shift condition. A given shift 
might eliminate the discontinuities in some components of the spatial metric, but it can 
not eliminate the shocks completely: at least some of the dynamical quantities will remain 
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non-smooth for all possible shift choices. 

Although in this paper I have concentrated in the BM hyperbolic formalism, the mathe- 
matical tools developed can easily be applied to any other hyperbolic formalism of General 
Relativity. One should expect the phenomena of coordinate shocks to also arise in any such 
formalism. In fact, since all formalisms must have the same physical solutions, the results 
of this paper imply that in any formalism the use of a harmonic slicing will generate shocks 
for at least some initial conditions. 

Clearly, the search for gauge conditions and/or restrictions on the initial data that can 
prevent the development of coordinate shocks is a problem that must be addressed if hy- 
perbolic formalisms are to become an important tool in the study of both theoretical and 
numerical relativity. 
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FIG. 1. Two-dimensional fiat spacetime. Initial values of the dynamical variables. 
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FIG. 3. Two-dimensional fiat spacetime. Values of the variables at t = 75 in the case when 
/ = 1.69. 
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FIG. 4. Two-dimensional fiat spacetime. Values of the variables at t = 75 in the case when 
/ = 0.49. 
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FIG. 6. Spherically symmetric flat spacetime. Initial values of the dynamical variables. 
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FIG. 7. Spherically symmetric flat spacetime. Values of the variables at t = 70 in the case 
when / = 1.69. 
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FIG. 8. Spherically symmetric flat spacetime. Values of the variables at t = 70 in the case 
when / = 0.49. 
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FIG. 9. Spherically symmetric flat spacetime. Values of the variables at t = 70 in the case 
when / = 1. 
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FIG. 10. Spherically symmetric black hole spacetime. Initial values of the dynamical variables. 
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FIG. 11. Spherically symmetric black hole spacetime. Values of the variables at t = 15 in the 
case when / = 1. 
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